{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Multi-Qubit Noisy Simulator\n",
    "*Copyright (c) 2021 Institute for Quantum Computing, Baidu Inc. All Rights Reserved.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Outline\n",
    "\n",
    "This tutorial will introduce how to use multi-qubit simulator at the pulse level. The outline is as follows:\n",
    "- Introduction\n",
    "- Preparation\n",
    "- Use multi-qubit noisy simulator at the gate level\n",
    "- Use multi-qubit noisy simulator at the pulse level\n",
    "    - Modeling the system \n",
    "    - Rabi oscillation\n",
    "    - Cross-Resonance effect\n",
    "    - ZZ crosstalk characterization through Ramsey experiment\n",
    "- Summary"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "\n",
    "Simulating time evolution of the qubits at the pulse level gives us more insight into the physics of quantum gates and the effects of noise. For superconducting quantum circuits, the transmon qubits are controlled by applying microwave pulses and magnetic flux. However, the performance of quantum gates is often suppressed by various factors: the decoherence of the qubit due to its interaction with the environment, the unwanted crosstalk effect, and leakage into the higher levels of the transmon. \n",
    "\n",
    "The multi-qubit noisy simulator provided by Quanlse allows us to simulate quantum operations on a noisy quantum device consisting of multiple transmon qubits to understand the physics behind quantum computing better. Several main types of noise are included in our noisy simulator: decoherence noise, amplitude noise, and crosstalk noise. We will focus on several common applications in superconducting quantum computing based on this noisy simulator: Rabi oscillation, Cross-Resonance effect, and characterizing ZZ crosstalk through a Ramsey experiment."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preparation\n",
    "\n",
    "After you have successfully installed Quanlse, you could run the Quanlse program below following this tutorial. To run this particular tutorial, you would need to import the following packages from Quanlse and other commonly-used Python libraries:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "is_executing": true
    }
   },
   "outputs": [],
   "source": [
    "from Quanlse.remoteOptimizer import remoteOptimize1Qubit as optimize1q\n",
    "from Quanlse.remoteSimulator import remoteSimulatorRunHamiltonian as runHamiltonian\n",
    "from Quanlse.Superconduct.Simulator import PulseModel\n",
    "from Quanlse.Superconduct.Simulator.PulseSim3Q import pulseSim3Q\n",
    "\n",
    "from Quanlse.QWaveform import QJob, QJobList\n",
    "from Quanlse.QOperator import driveX, driveY, a, sigmaZ, number, driveZ\n",
    "from Quanlse.QWaveform import square, gaussian\n",
    "from Quanlse.Utils.Functions import basis, tensor, expect, dagger, partialTrace, project, computationalBasisList, population\n",
    "from Quanlse.Utils.Bloch import rho2Coordinate, plotBloch\n",
    "from Quanlse.Utils.Plot import plotBarGraph\n",
    "from Quanlse.QOperation.FixedGate import H, X, Y, CR, Z, CNOT\n",
    "from Quanlse.QOperation.RotationGate import RZ\n",
    "from Quanlse.Superconduct.SchedulerSupport import SchedulerSuperconduct\n",
    "from Quanlse.Utils.Infidelity import unitaryInfidelity\n",
    "from Quanlse.Superconduct.SchedulerSupport.PipelineCenterAligned import centerAligned\n",
    "\n",
    "from math import pi\n",
    "from scipy.optimize import curve_fit\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To use Quanlse Cloud Service, we need to acquire a token to get access to the cloud. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "is_executing": true
    }
   },
   "outputs": [],
   "source": [
    "# Import Define class and set the token for cloud service\n",
    "# Please visit http://quantum-hub.baidu.com\n",
    "from Quanlse import Define\n",
    "Define.hubToken = \"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Use multi-qubit noisy simulator at the gate level\n",
    "\n",
    "We can also run the simulator at the gate level. In this section, we use a predefined `PulseModel()` instance with the default configuration. To create a 3-qubit physics model, we first instantiate the `PulseModel()` object by calling `pulseSim3Q()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = pulseSim3Q(frameMode='lab', dt=0.01)\n",
    "model.savePulse = False\n",
    "model.pipeline.addPipelineJob(centerAligned)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To define a circuit of creating GHZ state by Quanlse scheduler, we add the gates to the model by `gate(model.Q[index])`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Hadamard gate \n",
    "H(model.Q[0])\n",
    "\n",
    "# CNOT\n",
    "CNOT(model.Q[0], model.Q[1])\n",
    "CNOT(model.Q[1], model.Q[2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The pulse sequence of the quantum circuit defined above is generated by calling method `model.schedule`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scheJob = model.schedule()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the initial state $|000\\rangle$ and run the simulation. Then, plot the probability distribution of different outcome."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run the simulation\n",
    "res = model.simulate(job=scheJob)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the result\n",
    "popList = [abs(item ** 2) for item in res[0]['state'].T[0]]\n",
    "basisList = computationalBasisList(3, 3)\n",
    "plotBarGraph(basisList, popList, \"Result\", \"Outcome\", \"Population\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As it shows, the measurement result included unexpected values due to the noise. The simulation under decoherence can be done by setting the parameter `isOpen=True` in module `runHamiltonian()`, which takes some time to run, and we will obtain the density matrix after the simulation. For more details about how decoherence noise affects superconducting quantum computing, please refer to [Single-Qubit Noisy Simulator](https://quanlse.baidu.com/#/doc/tutorial-single-qubit-noisy-simulator)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Use multi-qubit noisy simulator at the pulse level\n",
    "\n",
    "Multi-qubit noisy simulator supports the quantum control simulation at the pulse level - get the system's final state by defining the waveform of the pulse and other related parameters. It allows us to simulate the quantum control of the superconducting hardware at the lower level. Here, we simulate some of the common operations used in real experiments.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Modeling the system\n",
    "\n",
    "Usually, we use Duffing oscillator to describe the physics model of the superconducting circuits. In lab frame, the system Hamiltonian of a three-qubit system with coupling between qubits $q_0$ and $q_1$, $q_1$ and $q_2$ reads:\n",
    "\n",
    "$$\n",
    "\\hat{H} = \\sum_{i=0}^2 \\omega_i \\hat{a}^\\dagger_i \\hat{a}_i + \\sum_{i=0}^2 \\frac{\\alpha_i}{2} \\hat{a}^\\dagger_i \\hat{a}^\\dagger_i \\hat{a}_i \\hat{a}_i + g_{01} (\\hat{a}^\\dagger_0 \\hat{a}_1 + \\hat{a}_0 \\hat{a}^\\dagger_1) +  g_{12} (\\hat{a}^\\dagger_1 \\hat{a}_2 + \\hat{a}_1 \\hat{a}^\\dagger_2),\n",
    "$$\n",
    "\n",
    "where $\\omega_i$ and $\\alpha_i$ are the qubit frequency and anharmonicity of qubit $q_i$ respectively; $g_{i, j}$ is the coupling strength between qubit $q_i$ and qubit $q_j$; $a_i$, $a^\\dagger_i$ denote the annihilation operator and creation operator of qubit $q_i$.\n",
    "\n",
    "In this tutorial, we use a three-qubit system as an example. We first define the parameters of the hardware."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "is_executing": true
    }
   },
   "outputs": [],
   "source": [
    "qubitNum = 3  # The number of qubits\n",
    "level = 3  # The energy level for each qubit\n",
    "\n",
    "anharm = -0.33 * 2 * pi  # The anharmonicity of the qubit, in 2 pi GHz\n",
    "wq0 = 4.914 * 2 * pi  # The frequency for qubit 0, in 2 pi GHz \n",
    "wq1 = 5.100 * 2 * pi  # The frequency for qubit 1, in 2 pi GHz\n",
    "wq2 = 5.200 * 2 * pi  # The frequency for qubit 2, in 2 pi GHz\n",
    "g01 = 0.0038 * 2 * pi  # The coupling strength of the interaction between qubit 0 and qubit 1, in 2 pi GHz\n",
    "g12 = 0.0020 * 2 * pi  # The coupling strength of the interaction between qubit 1 and qubit 2, in 2 pi GHz\n",
    "\n",
    "dt = 1.  # The sampling time of AWG\n",
    "\n",
    "# T1 relaxation time for qubit 0, qubit 1, and qubit 2, in nanoseconds\n",
    "t01 = 1000  \n",
    "t11 = 1120\n",
    "t21 = 1300\n",
    "\n",
    "# T2 dephasing time for qubit 0, qubit 1, and qubit 2, in nanoseconds\n",
    "t02 = 500\n",
    "t12 = 450\n",
    "t22 = 600\n",
    "\n",
    "# The random amplitude distortion\n",
    "ampNoise = 0.02"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The physics model is created by instantiating an object of class `PulseModel`. The types of noise include $T_1$-relaxation noise, $T_2$-dephasing, and distortion of amplitudes.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "is_executing": true
    }
   },
   "outputs": [],
   "source": [
    "qubitFreq = {0: wq0, 1: wq1, 2: wq2}  # Qubit frequency for each qubit\n",
    "qubitAnharm = {0: anharm, 1: anharm, 2: anharm}  # Qubit anharmonicity for each qubit\n",
    "qubitT1 = {0: t01, 1: t11, 2: t21}  # Relaxation time \n",
    "qubitT2 = {0: t02, 1: t12, 2: t22}  # Dephasing time\n",
    "couplingMap = {(0, 1): g01, (1, 2): g12}  # Coupling map\n",
    "\n",
    "# Create an instant of PulseModel\n",
    "model = PulseModel(subSysNum=qubitNum,\n",
    "                   sysLevel=level,\n",
    "                   qubitFreq=qubitFreq,\n",
    "                   qubitAnharm=qubitAnharm,\n",
    "                   couplingMap=couplingMap,\n",
    "                   T1=qubitT1,\n",
    "                   T2=qubitT2,\n",
    "                   dt=dt,\n",
    "                   ampSigma=ampNoise)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have constructed a noisy simulator including three superconducting qubits with three types of noises. The next step is to create a `QHamiltonian` object by calling method `createQHamiltonian()`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "is_executing": true
    }
   },
   "outputs": [],
   "source": [
    "ham = model.createQHamiltonian()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cross-Resonance effect \n",
    "\n",
    "The all-microwave control is one of the strategies to realize quantum control on superconducting circuits. In this strategy, two-qubit operations harness the cross-resonance effect of two weakly-coupled qubits. This is done by driving the control qubit with the frequency of the weakly-coupled target qubit. Ideally, the desired $\\hat{\\sigma}_z \\otimes \\hat{\\sigma}_x$ interaction between the control and target qubit is dominating interaction \\[1\\]. For more details about CR gate, please refer to [Cross-Resonance Gate](https://quanlse.baidu.com/#/doc/tutorial-cr).\n",
    "\n",
    "In our simulation, we again drive qubit $q_0$ (control qubit) by various amplitudes (with the drive frequency of the qubit $q_1$). This can be done by `addWaveRot(index, waves, detuning)` where `index` is the index of qubit acted upon; `waves` is the waveform of the pulse; and `detuning` $\\Delta$ is the frequency difference ($\\Delta = \\omega_q - \\omega_d$, where $\\omega_q$ is the qubit frequency and $\\omega_d$ the drive frequency). \n",
    "\n",
    "Here, we vary the amplitudes of the pulse and record the population of $|1\\rangle$ for each qubit. In this example, $q_1$ is the control qubit driven by the pulse with the frequency of target qubit $q_0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dCoef = 0.03 * (2 * pi)  # The drive strength of the pulse\n",
    "ampCR = np.linspace(0, 0.5, 40)  # The amplitudes in arbitrary unit \n",
    "amps = ampCR * dCoef  \n",
    "detuning = wq1 - wq0  # The detuning of the pulse\n",
    "\n",
    "# jobList = QJobList(subSysNum=qubitNum, sysLevel=level, dt=dt, title='cr')\n",
    "jobList = ham.createJobList()\n",
    "\n",
    "# Fix the gate time\n",
    "tg = 950\n",
    "\n",
    "# Append each job to the jobList\n",
    "for amp in amps:\n",
    "    job = QJob(subSysNum=qubitNum, sysLevel=level, dt=dt)\n",
    "    job = ham.createJob()\n",
    "    job.addWaveRot(1, waves=square(tg, amp), t0=0., detuning=detuning)  # Apply pulse at qubit 1\n",
    "    job = model.getSimJob(job)\n",
    "    jobList.addJob(jobs=job)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the simulation with initial state $|\\psi\\rangle = |010\\rangle$, and the control qubit $q_1$ starts at the excited state."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "is_executing": true
    }
   },
   "outputs": [],
   "source": [
    "# Define the initial state of |010>\n",
    "psi0 = tensor(basis(level, 0), basis(level, 1), basis(level, 0))  \n",
    "\n",
    "# Run the simulation\n",
    "result = runHamiltonian(ham=ham, state0=psi0, jobList=jobList)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the projector of the first excited state for qubit $q_i$ and initialize the list of expected values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "is_executing": true
    }
   },
   "outputs": [],
   "source": [
    "prj01 = tensor(basis(3, 1) @ dagger(basis(3,1)), np.identity(level), np.identity(level))  # The projector of qubit 0\n",
    "prj11 = tensor(np.identity(level), basis(3, 1) @ dagger(basis(3,1)), np.identity(level))  # The projector of qubit 1\n",
    "prj21 = tensor(np.identity(level), np.identity(level), basis(3, 1) @ dagger(basis(3,1)))  # The projector of qubit 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute the expected values of the projector of each qubit, and plot them with respect to the different amplitudes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "is_executing": true
    }
   },
   "outputs": [],
   "source": [
    "# Initialize the list of expected values\n",
    "num0List = []\n",
    "num1List = []\n",
    "num2List = []\n",
    "\n",
    "for res in result.result:\n",
    "    state = res['state']  # The final state of each job\n",
    "    num0Expect = expect(prj01, state)  # Compute the expected values of the projector |1><1|\n",
    "    num1Expect = expect(prj11, state)\n",
    "    num2Expect = expect(prj21, state)\n",
    "    num0List.append(num0Expect)\n",
    "    num1List.append(num1Expect)\n",
    "    num2List.append(num2Expect)\n",
    "\n",
    "plt.figure(figsize=(8, 6))\n",
    "# Plot the figure of CR effect\n",
    "plt.plot(ampCR, num0List, label='qubit0')\n",
    "plt.plot(ampCR, num1List, label='qubit1')\n",
    "plt.plot(ampCR, num2List, label='qubit2')\n",
    "\n",
    "plt.xlabel('Amplitudes (a.u.)')\n",
    "plt.ylabel(r'Population of $|1\\rangle$')\n",
    "plt.title('Cross-Resonance effect')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As it shows, the projector $|1\\rangle \\langle 1|$'s expected value of target qubit $q_0$ changes over the amplitude, while the control qubit $q_1$ is in the excited state when the amplitude is relatively small. Meanwhile, qubit $q_2$ is always in the ground state. It can also be seen that the increasing amplitude inevitably affects qubit $q_1$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ZZ crosstalk characterization by a Ramsey experiment\n",
    "\n",
    "ZZ crosstalk is the major source of unwanted interaction between coupled qubits. It arises from the existence of states of higher energy levels. The effective Hamiltonian of two coupled qubits (directly or indirectly) in the two-qubit subspace is \\[2\\]:\n",
    "\n",
    "$$\n",
    "\\hat{H}_{\\rm eff} = \\omega_{0}\\frac{\\hat{\\sigma}_{z}^0 \\otimes I_1}{2} + \\omega_{1}\\frac{I_0\\otimes\\hat{\\sigma}_{z}^1}{2} + \\xi \\frac{\\hat{\\sigma}_{z}^0 \\otimes \\hat{\\sigma}_{z}^1}{2},\n",
    "$$\n",
    "\n",
    "Where $\\omega_0$, $\\omega_1$ are the qubit frequencies and $\\xi$ is the strength of ZZ crosstalk. $\\xi$ is defined as the different of transition frequencies between $|11\\rangle \\leftrightarrow |10\\rangle$ and $|01\\rangle \\leftrightarrow |00\\rangle$:\n",
    "\n",
    "$$\n",
    "\\xi = \\left(E_{11} - E_{10}\\right) - \\left(E_{01} - E_{00}\\right),\n",
    "$$\n",
    "\n",
    "where $E_{ij}$ is the energy level of state $|ij\\rangle$. We can actually detect and measure this frequency shift-induced crosstalk by Ramsey experiment. This can be done by applying two Hadamard gates with an idle time apart \\[3\\]. \n",
    "\n",
    "To better illustrate the effect of ZZ crosstalk, we define a new 3-qubit model with stronger coupling strengths (6 ~ 40 MHz)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dt = 0.2  # The sampling time\n",
    "level = 3  # The system level\n",
    "qubitNum = 3  # The number of qubits\n",
    "\n",
    "g01 = 0.0377 * (2 * pi)\n",
    "g12 = 0.0060 * (2 * pi)\n",
    "\n",
    "# Coupling map\n",
    "couplingMap = {\n",
    "    (0, 1): g01,\n",
    "    (1, 2): g12\n",
    "}\n",
    "\n",
    "# Qubits frequency anharmonicity\n",
    "anharm = - 0.33 * (2 * pi)\n",
    "qubitAnharm = {0: anharm, 1: anharm, 2: anharm}  # The anharmonicities for each qubit\n",
    "\n",
    "# Qubit Frequency\n",
    "qubitFreq = {\n",
    "            0: 5.5904 * (2 * pi),\n",
    "            1: 4.7354 * (2 * pi),\n",
    "            2: 4.8524 * (2 * pi)\n",
    "            }"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create the physics model by class `PulseModel()`, and create Hamiltonian `ham` by the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = PulseModel(subSysNum=qubitNum, sysLevel=level, couplingMap=couplingMap,\n",
    "                    qubitFreq=qubitFreq, dt=dt, qubitAnharm=qubitAnharm)\n",
    "    \n",
    "ham = model.createQHamiltonian()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate the pulses of the gates $H$ and $X$ on different qubits."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define function to generate the QJob for gate of specified qubit \n",
    "def generateGate(gate, index):\n",
    "    job1q, _ = optimize1q(ham=ham.subSystem(index), uGoal=gate.getMatrix(), targetInfid=1e-5)\n",
    "    job3q = QJob(subSysNum=qubitNum, sysLevel=level, dt=dt)\n",
    "    waves = job1q.waves\n",
    "    ops = job1q.ctrlOperators\n",
    "    for key, op in ops.items():\n",
    "        job3q.addWave(operators=op, onSubSys=index, waves=waves[key])\n",
    "     \n",
    "    return job3q\n",
    "\n",
    "# Generate the gates needed\n",
    "h0 = generateGate(H, 0)  # H gate on qubit 0 \n",
    "h1 = generateGate(H, 1)  # H gate on qubit 1\n",
    "x1 = generateGate(X, 1)  # X gate on qubit 1\n",
    "x2 = generateGate(X, 2)  # X gate on qubit 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "maxTime = 500  # The delayed time in Ramsey experiment, in nanosecond.\n",
    "freq = 3 / maxTime  # Detuning. \n",
    "\n",
    "# Generate job for delayed time \n",
    "def generateIdle(tg, index):\n",
    "    jobIdle = QJob(subSysNum=qubitNum, sysLevel=level, dt=dt)\n",
    "    jobIdle.appendWave(operators=driveZ, onSubSys=index, waves=square(tg, 2 * pi * freq))\n",
    "    return jobIdle"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define two different `jobList` objects - one begins with $|00\\rangle$ and the other $|01\\rangle$ by applying a $X$ gate on qubit $q_1$. Then perform Ramsey experiment on qubit $q_0$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# jobList with initial state |00>\n",
    "jobListGrd = QJobList(subSysNum=qubitNum, sysLevel=level, dt=dt)\n",
    "\n",
    "# jobList with initial state |01> (by applying X gate) \n",
    "jobListExd = QJobList(subSysNum=qubitNum, sysLevel=level, dt=dt)\n",
    "\n",
    "# Define the delayed time\n",
    "tgList = np.linspace(0, maxTime, 50)\n",
    "\n",
    "# Define jobList with initial state |00>\n",
    "for tg in tgList:\n",
    "    job = QJob(subSysNum=qubitNum, sysLevel=level, dt=dt)\n",
    "    job += h0\n",
    "    job += generateIdle(tg, 0)\n",
    "    job += h0\n",
    "    jobListGrd.addJob(job)\n",
    "\n",
    "# Define jobList with initial state |01>\n",
    "for tg in tgList:\n",
    "    job = QJob(subSysNum=qubitNum, sysLevel=level, dt=dt)\n",
    "    job += x1\n",
    "    job += h0\n",
    "    job += generateIdle(tg, 0)\n",
    "    job += h0\n",
    "    jobListExd.addJob(job)\n",
    "\n",
    "# Run the simulation\n",
    "stateInit = tensor(basis(level, 0), basis(level, 0), basis(level, 0))\n",
    "resultGrd = runHamiltonian(ham, state0=stateInit, jobList=jobListGrd)\n",
    "resultExd = runHamiltonian(ham, state0=stateInit, jobList=jobListExd)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the population of excited state of qubit $q_0$ versus delayed time. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "num0List = []\n",
    "num1List = []\n",
    "\n",
    "# projector |1><1| of qubit 0\n",
    "prj1 = tensor(basis(level, 1) @ dagger(basis(level, 1)), np.identity(9))\n",
    "\n",
    "# append the result to the list\n",
    "for res0, res1 in zip(resultGrd, resultExd):\n",
    "    psi0, psi1 = res0['state'], res1['state']\n",
    "    rho0, rho1 = psi0 @ dagger(psi0), psi1 @ dagger(psi1)\n",
    "    num0List.append(expect(prj1, rho0))\n",
    "    num1List.append(expect(prj1, rho1))\n",
    "\n",
    "# plot the result\n",
    "plt.figure(figsize=(8, 6))\n",
    "plt.plot(tgList, num0List, '.b', label=r'$|00\\rangle$')\n",
    "plt.plot(tgList, num1List, '.r', label=r'$|01\\rangle$')\n",
    "plt.xlabel('Delayed time (ns)')\n",
    "plt.ylabel('Population of excited state of qubit 0')\n",
    "plt.xlabel\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The strength of ZZ crosstalk $\\xi$ can be estimated by computing the frequency difference in the Ramsey oscillation. Therefore, we use the cosine function to fit the result acquired by simulation to compute the frequencies $f_1$, $f_2$. The strength is given by $\\xi / \\left( 2\\pi \\right) = |f_1 - f_2|$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the fitting curve\n",
    "def fit(x, omega, theta):\n",
    "    return - 0.5 * np.cos(omega * x + theta) + 0.5\n",
    "\n",
    "# Fit the curve\n",
    "para1Fit, _ = curve_fit(fit, tgList, num0List, [2.1 * pi * freq, 0])\n",
    "para2Fit, _ = curve_fit(fit, tgList, num1List, [2 * pi * freq, 0])\n",
    "step = 0.01\n",
    "y1Fit = [fit(x, para1Fit[0], para1Fit[1]) for x in np.arange(tgList[0], tgList[-1], step)]\n",
    "y2Fit = [fit(x, para2Fit[0], para2Fit[1]) for x in np.arange(tgList[0], tgList[-1], step)]\n",
    "\n",
    "# Plot the curve\n",
    "plt.figure(figsize=(8, 6))\n",
    "plt.plot(np.arange(tgList[0], tgList[-1], step), y1Fit)\n",
    "plt.plot(np.arange(tgList[0], tgList[-1], step), y2Fit)\n",
    "plt.plot(tgList, num0List, '.b', label=r'$|00\\rangle$')\n",
    "plt.plot(tgList, num1List, '.r', label=r'$|01\\rangle$')\n",
    "plt.xlabel('Delayed time (ns)')\n",
    "plt.ylabel('Population of excited state of qubit 0')\n",
    "plt.title('Ramsey on Q0')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate the crosstalk strength\n",
    "xiEst = abs(para1Fit[0] - para2Fit[0]) \n",
    "print(f'Coupling strength: {g01 * 1e3 / (2 * pi)} MHz')\n",
    "print(f'ZZ crosstalk strength: {xiEst * 1e3 / (2 * pi)} MHz')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Due to the strong coupling strength between qubits $𝑞_0$ and $𝑞_1$, it can be observed that the frequency difference is relatively large, that is, the 𝑍𝑍 crosstalk is relatively large."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can repeat the same experiment to calculate ZZ crosstalk strength $\\xi$ between qubit $q_1$ and qubit $q_2$ with smaller coupling strength."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# jobList with initial state |00>\n",
    "jobListGrd = QJobList(subSysNum=qubitNum, sysLevel=level, dt=dt)\n",
    "\n",
    "# jobList with initial state |01> (by applying X gate)\n",
    "jobListExd = QJobList(subSysNum=qubitNum, sysLevel=level, dt=dt)\n",
    "\n",
    "# Define the delayed time\n",
    "tgList = np.linspace(0, maxTime, 50)\n",
    "\n",
    "# Define jobList with initial state |00>\n",
    "for tg in tgList:\n",
    "    job = QJob(subSysNum=qubitNum, sysLevel=level, dt=dt)\n",
    "    job += h1\n",
    "    job += generateIdle(tg, 1)\n",
    "    job += h1\n",
    "    jobListGrd.addJob(job)\n",
    "\n",
    "# Define jobList with initial state |01>    \n",
    "for tg in tgList:\n",
    "    job = QJob(subSysNum=qubitNum, sysLevel=level, dt=dt)\n",
    "    job += x2\n",
    "    job += h1\n",
    "    job += generateIdle(tg, 1)\n",
    "    job += h1\n",
    "    jobListExd.addJob(job)\n",
    "\n",
    "# Run the simulation    \n",
    "stateInit = tensor(basis(level, 0), basis(level, 0), basis(level, 0))\n",
    "resultGrd = runHamiltonian(ham, state0=stateInit, jobList=jobListGrd)\n",
    "resultExd = runHamiltonian(ham, state0=stateInit, jobList=jobListExd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "num0List = []\n",
    "num1List = []\n",
    "\n",
    "# projector |1><1| of qubit 1\n",
    "prj1 = tensor(np.identity(3), basis(level, 1) @ dagger(basis(level, 1)), np.identity(3))\n",
    "\n",
    "# append the result to the list\n",
    "for res0, res1 in zip(resultGrd, resultExd):\n",
    "    psi0, psi1 = res0['state'], res1['state']\n",
    "    rho0, rho1 = psi0 @ dagger(psi0), psi1 @ dagger(psi1)\n",
    "    num0List.append(expect(prj1, rho0))\n",
    "    num1List.append(expect(prj1, rho1))\n",
    "\n",
    "\n",
    "# plot the result\n",
    "plt.figure(figsize=(8, 6))\n",
    "plt.plot(tgList, num0List, '.b', label=r'$|00\\rangle$')\n",
    "plt.plot(tgList, num1List, '.r', label=r'$|01\\rangle$')\n",
    "plt.xlabel('Delayed time (ns)')\n",
    "plt.ylabel('Population of excited state of qubit 1')\n",
    "plt.xlabel\n",
    "plt.title('Ramsey on Q1')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Fit the curve\n",
    "para1Fit, _ = curve_fit(fit, tgList, num0List, [2 * pi * freq / 1.2, 0.])\n",
    "para2Fit, _ = curve_fit(fit, tgList, num1List, [2 * pi * freq / 1.2, 0.])\n",
    "step = 0.01\n",
    "y1Fit = [fit(x, para1Fit[0], para1Fit[1]) for x in np.arange(tgList[0], tgList[-1], step)]\n",
    "y2Fit = [fit(x, para2Fit[0], para2Fit[1]) for x in np.arange(tgList[0], tgList[-1], step)]\n",
    "\n",
    "# Plot the curve\n",
    "plt.figure(figsize=(8, 6))\n",
    "plt.plot(np.arange(tgList[0], tgList[-1], step), y1Fit)\n",
    "plt.plot(np.arange(tgList[0], tgList[-1], step), y2Fit)\n",
    "plt.plot(tgList, num0List, '.b', label=r'$|00\\rangle$')\n",
    "plt.plot(tgList, num1List, '.r', label=r'$|01\\rangle$')\n",
    "plt.xlabel('Delayed time (ns)')\n",
    "plt.ylabel('Population of excited state of qubit 1')\n",
    "plt.xlabel\n",
    "plt.title('Ramsey on Q1')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate the crosstalk strength\n",
    "xiEst = abs(para1Fit[0] - para2Fit[0]) \n",
    "print(f'Coupling strength: {g12 * 1e3 / (2 * pi)} MHz')\n",
    "print(f'ZZ crosstalk strength: {xiEst * 1e3 / (2 * pi)} MHz')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Due to the weaker coupling strength, the relatively small qubit frequency shift of $q_1$ indicates the weak ZZ crosstalk between qubit $q_1$ and $q_2$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summary\n",
    "\n",
    "After reading this tutorial on multi-qubit noisy simulator, the users could follow this link [tutorial-multi-qubit-noisy-simulator.ipynb](https://github.com/baidu/Quanlse/blob/main/Tutorial/EN/tutorial-multi-qubit-noisy-simulator.ipynb) to the GitHub page of this Jupyter Notebook document and run this program for themselves. The users are encouraged to explore other advanced research which is different from this tutorial."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\\[1\\] [Malekakhlagh, Moein, Easwar Magesan, and David C. McKay. \"First-principles analysis of cross-resonance gate operation.\" *Physical Review A* 102.4 (2020): 042605.](https://journals.aps.org/pra/abstract/10.1103/PhysRevA.102.042605)\n",
    "\n",
    "\\[2\\] [Magesan, Easwar, and Jay M. Gambetta. \"Effective Hamiltonian models of the cross-resonance gate.\" *Physical Review A* 101.5 (2020): 052308.](https://journals.aps.org/pra/abstract/10.1103/PhysRevA.101.052308)\n",
    "\n",
    "\\[3\\] [Ku, Jaseung, et al. \"Suppression of Unwanted ZZ Interactions in a Hybrid Two-Qubit System.\" *Physical review letters* 125.20 (2020): 200504.](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.125.200504)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
